﻿
using System;
using PVTLibrary.Packages;

namespace PVTLibrary.FlashAlgorithms
{
    [Serializable] 
    public class BostonBrittInsideOut : FlashAlgorithm
    {


        private int i, j, k, n, ecount;
        double etol = 0.000001;
        double itol = 0.000001;
        private int maxit_i = 100;
        int maxit_e  = 100;
        string[] Vn;
        double[] Vx, Vy, Vp, ui, uic, pi, Ki, fi, Vt, Vpc, VTc, Vw;
        double L, Lf, V, Vf, R, Rant, Tant, Pant, T, T_, Tf, P, P_, Pf, T0, P0, A, B, C, D, E, F, Ac, Bc, Cc, Dc, Ec, Fc;
        private double Kb, Kb0, Kb_;
        private double DHv, DHl, DHv1, DHv2, DHl1, DHl2, Hv0, Hvid, Hlid, Hf, DHlsp, DHvsp;
        private double DSv, DSl, DSv1, DSv2, DSl1, DSl2, Sv0, Svid, Slid, Sf, DSlsp, DSvsp;
        private double Pb, Pd, Pmin, Pmax, Px, soma_x, soma_y, Tmin, Tmax;
        private PropertyPackage proppack;
        private object tmpdx, refx, currx;


        /// <summary>
        /// 
        /// </summary>
        /// <param name="Vz"></param>
        /// <param name="P"></param>
        /// <param name="T"></param>
        /// <param name="PP"></param>
        /// <param name="ReuseKI"></param>
        /// <param name="PrevKi"></param>
        public override object[] Flash_PT(double[] Vz , double P , double T , PropertyPackage PP , bool ReuseKI  = false, double[] PrevKi =null)
        {
            DateTime d1, d2;
            TimeSpan dt;

            d1 = DateTime.Now;

            etol = (double) PP.Parameters("PP_PTFELT");
            maxit_e = (int) PP.Parameters("PP_PTFMEI");
            itol = (double) PP.Parameters("PP_PTFILT");
            maxit_i = (int) PP.Parameters("PP_PTFMII");


            proppack = PP;

            n = Vz.Length;
            Vn = new string[n];
            Vx = new double[n];
            Vy = new double[n];
            Vp = new double[n];

            ui = new double[n];
            uic = new double[n];
            pi = new double[n];
            Ki = new double[n];
            fi = new double[n];


            Vn = PP.RET_VNAMES();
            fi = Vz.Clone();

            //--------------------------------------
            // STEP 1 - Assume u, A, B, C, D, E, F 
            //--------------------------------------


            //----------------------------------------
            // STEP 1.1 - Estimate K, Vx, Vy, V and L 
            //----------------------------------------

            // Calculate Ki`s


            if(!ReuseKI)
            {
                i = 0;
                do
                {
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T);
                    Ki[i] = Vp[i]/P;
                    i += 1;
                } while (i < n + 1);
            }
            else
            {
                for( i = 0 ;i< n;i++)
                {
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T);
                    Ki[i] = PrevKi[i];
                }
            }

            ////Estimate V

            if( T > Common.Max(proppack.RET_VTC, Vz))
            {
                Vy = Vz;
                V  = 1;
                L  = 0;
                goto out;
            }

            i  = 0;
            Px = 0;
            do
            {
                Px = Px + (Vz[i]/Vp[i]);
                i = i + 1;
            } while (i < n + 1);
            Px = 1/Px;
            Pmin = Px;
            i = 0;
            Px = 0;

            do
            {
                Px = Px + Vz[i]*Vp[i];
                i = i + 1;
            } while (i < n + 1);

            Pmax = Px;
            Pb = Pmax;
            Pd = Pmin;

            if( Math.Abs(Pb - Pd) / Pb < 0.0000001 )
            {
                //one comp only
                if( Px <= P)
                {
                    L = 1;
                    V = 0;
                    Vx = Vz;
                    goto out;
                }
                else
                {
                    L = 0;
                    V = 1;
                    Vy = Vz;
                    goto out;
                }
            }
            else if (P <= Pd)
            {
                //vapor only
                L = 0.02;
                V = 0.98;
            }
            else if (P >= Pb)
            {
                //liquid only
                L = 0.02;
                V = 0.98;
            }
            else
            {
                //VLE
                V = 1 - (P - Pd)/(Pb - Pd);
                L = 1 - V;
            }

            if(n == 0)
            {
                if( Vp(0) <= P)
                {
                    L = 1;
                    V = 0;
                    Vx = Vz;
                    goto out;
                }
                else
                {
                    L = 0;
                    V = 1;
                    Vy = Vz;
                    goto out;
                }
            }

            i = 0;
            do
            {
                if (Vz[i] != 0)
                {
                    Vy[i] = Vz[i]*Ki[i]/((Ki[i] - 1)*V + 1);
                    Vx[i] = Vy[i]/Ki[i];
                }
                else
                {
                    Vy[i] = 0;
                    Vx[i] = 0;
                }
                i += 1
            } while (i < n + 1);

            i = 0;
            soma_x = 0;
            soma_y = 0;
            do
            {
                soma_x = soma_x + Vx[i];
                soma_y = soma_y + Vy[i];
                i = i + 1;
            } 
            while (i < n + 1);

            i = 0;
            do
            {
                Vx[i] = Vx[i]/soma_x;
                Vy[i] = Vy[i]/soma_y;
                i = i + 1;
            } 
            while (i < n + 1);

            Kb = 1;//CalcKbj1(Ki)
            Kb0 = 1;//Kb

            for( i = 0 ;i< n;i++)
            {
                ui[i] = Math.Log(Ki[i]);
            }

            double[] fx, x, dx, xbr, fbr;
            double[,] dfdx;
            Brent bo = new Brent();
            BrentMinimize bo2 = new BrentMinimize();

            ecount = 0;

            do
            {
                //--------------------------------------------------------------
                // STEPS 2, 3, 4, 5, 6, 7 and 8 - Calculate R and Energy Balance
                //--------------------------------------------------------------
                Rant = R;
                R = Kb * V / (Kb * V + Kb0 * L);

                //Dim fr As Double

                //bo2.DefineFuncDelegate(AddressOf TPErrorFunc)
                //fr = bo2.brentoptimize(0, 1, 0.0000000001, R)

                double fr, dfr, R0, R1;
                int icount = 0;

                do
                {
                    R0 = R;
                    if (R > 0.999)
                    {
                        R1 = R - 0.001;
                        fr = TPErrorFunc(R0);
                        dfr = (fr - TPErrorFunc(R1))/0.001;
                    }
                    else
                    {
                        R1 = R + 0.001;
                        fr = TPErrorFunc(R0);
                        dfr = (fr - TPErrorFunc(R1))/-0.001;
                    }
                    R0 = R;

                    R = R - 0.3*fr/dfr;

                    if (R < 0.0)
                        R = 0.0;

                    if (R > 1.0)
                        R = 1.0;

                    icount += 1;
                } 
                while (Math.Abs(fr) > itol && icount < maxit_i + 1);

                if( icount > maxit_i)
                     R = Rant;

                if (Rant = 0.0 && R = 1.0)
                    R = 0.0;

                if (Rant = 1.0 && R = 0.0)
                    R = 1.0;

                TPErrorFunc(R);

                if( R <= 0)
                {
                    R = 0;
                    L = 1;
                    V = 0;
                    for(i = 0 ;i< n;i++)
                    {
                        pi[i] = fi[i]/(1 - R + Kb0*R*Math.Exp(ui[i]));
                    }
                    double sumpi = 0;
                    double sumeuipi = 0;
                    
                    for( i = 0;i< n;i++)
                    {
                        sumpi += pi[i];
                        sumeuipi += Math.Exp(ui[i])*pi[i];
                    }
                    
                    for( i = 0;i< n;i++)
                    {
                        Vx[i] = pi[i]/sumpi;
                        Vy[i] = Math.Exp(ui[i])*pi[i]/sumeuipi;
                    }
                    break;
                }

                if( R >= 1)
                {
                    R = 1;
                    L = 0;
                    V = 1;
                    for( i = 0;i< n;i++)
                    {
                        pi[i] = fi[i]/(1 - R + Kb0*R*Math.Exp(ui[i]));
                    }
                    double sumpi = 0;
                    Double sumeuipi = 0;
                    for( i = 0;i< n;i++)
                    {
                        sumpi += pi[i];
                        sumeuipi += Math.Exp(ui[i])*pi[i];
                    }
                    for( i = 0;i< n;i++)
                    {
                        Vx[i] = pi[i]/sumpi;
                        Vy[i] = Math.Exp(ui[i])*pi[i]/sumeuipi;
                    }
                    break;
                }


                //At this point, we have converged R for the simplified model. Proceed to step 9.

                //----------------------------------------------------------
                // STEP 9 - Rigorous model Enthalpy and K-values calculation
                //----------------------------------------------------------

                Ki = PP.DW_CalcKvalue(Vx, Vy, T, P);
                //Kb = CalcKbj1(Ki)

                for( i = 0;i< n;i++)
                {
                    uic[i] = Log(Ki[i])
                }

                //-------------------------------------------
                // STEP 10 - Update variables using Broyden
                //-------------------------------------------

                for( i = 0;i< n;i++)
                {
                    fx[i] = (ui[i] - uic[i]);
                    x(i) = ui[i];
                }

                if(ecount == 0)
                {
                    for( i = 0;i< n;i++)
                    {
                        for(j = 0 ;j<n;j++)
                        {
                            if (i == j)
                                dfdx[i][j] = 1;
                            else
                                dfdx[i][j] = 0;
                        }
                    }
                    Broyden.broydn(n, x, fx, dx, xbr, fbr, dfdx, 0);
                }
                else
                {
                    Broyden.broydn(n, x, fx, dx, xbr, fbr, dfdx, 1);
                }

                for( i = 0;i< n;i++)
                {
                    ui[i] = ui[i] + dx(i);
                }

                ecount += 1;

                If Double.IsNaN(V) 
                    Then Throw New Exception("Error calculating the vapor fraction.")
                
                If ecount > maxit_e 
                     Throw New Exception(DTL.App.GetLocalString("PropPack_FlashMaxIt2"))

                Console.WriteLine("PT Flash [IO]: Iteration #" & ecount & ", VF = " & V)


            } while ( AbsSum(fx) > etol && Math.Abs(R - Rant) > 0.0000000001))

            d2 = Date.Now;

            dt = d2 - d1;

            Console.WriteLine("PT Flash [IO]: Converged in " & ecount & " iterations. Time taken: " & dt.TotalMilliseconds & " ms. Error function value: " & AbsSum(fx))

out:        
            return new Object() {L, V, Vx, Vy, ecount, 0.0, Vx, 0.0#, PP.RET_NullVector}

            }






            public Overrides Function Flash_PH(ByVal Vz As Double(), ByVal P As Double, ByVal H As Double, ByVal Tref As Double, ByVal PP As PropertyPackages.PropertyPackage, Optional ByVal ReuseKI As Boolean = False, Optional ByVal PrevKi As Double() = Nothing) As Object

            Dim d1, d2 As Date, dt As TimeSpan

            d1 = Date.Now

            maxit_i = CInt(PP.Parameters("PP_PHFMII"))
            maxit_e = CInt(PP.Parameters("PP_PHFMEI"))
            itol = CDbl(PP.Parameters("PP_PHFILT"))
            etol = CDbl(PP.Parameters("PP_PHFELT"))

            n = UBound(Vz)

            proppack = PP
            Hf = H * PP.AUX_MMM(Vz)
            Pf = P

            ReDim Vn(n), Vx(n), Vy(n), Vp(n), ui(n), uic(n), pi(n), Ki(n), fi(n), Vpc(n), VTc(n), Vw(n)

            Vn = PP.RET_VNAMES()
            VTc = PP.RET_VTC
            fi = Vz.Clone

            Tmin = 0
            Tmax = 0
            If Tref = 0.0# Or Double.IsNaN(Tref) Then
                i = 0
                Tref = 0
                Do
                    Tref += 0.8 * Vz[i] * VTc(i)
                    Tmin += 0.1 * Vz[i] * VTc(i)
                    Tmax += 2.0 * Vz[i] * VTc(i)
                    i += 1
                Loop Until i = n + 1
            Else
                Tmin = Tref - 200
                Tmax = Tref + 200
            End If
            If Tmin < 100 Then Tmin = 100

            //--------------------------------------
            // STEP 1 - Assume u, A, B, C, D, E, F 
            //--------------------------------------

            T = Tref + 1
            T_ = Tref - 1
            T0 = Tref

            //----------------------------------------
            // STEP 1.1 - Estimate K, Vx, Vy, V and L 
            //----------------------------------------

            //Calculate Ki`s

            Vpc = PP.RET_VPC
            VTc = PP.RET_VTC
            Vw = PP.RET_VW

            If Not ReuseKI Then
                i = 0
                Do
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    Ki[i]= Vpc(i) / P * Math.Exp(5.373 * (1 + Vw(i)) * (1 - VTc(i) / T))
                    i += 1
                Loop Until i = n + 1
            Else
                for( i = 0;i< n;i++)
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    Ki[i]= PrevKi[i]
                Next
            End If

            //Estimate V

            i = 0
            Px = 0
            Do
                Px = Px + (Vz[i] / Vp[i])
                i = i + 1
            Loop Until i = n + 1
            Px = 1 / Px
            Pmin = Px
            i = 0
            Px = 0
            Do
                Px = Px + Vz[i] * Vp[i]
                i = i + 1
            Loop Until i = n + 1

            Pmax = Px
            Pb = Pmax
            Pd = Pmin

            If P <= Pd Then
                //vapor only
                L = 0
                V = 1
            ElseIf P >= Pb Then
                //liquid only
                L = 1
                V = 0
            Else
                //VLE
                V = 1 - (P - Pd) / (Pb - Pd)
                L = 1 - V
            End If

            If n = 0 Then
                If Vp(0) <= P Then
                    L = 1
                    V = 0
                Else
                    L = 0
                    V = 1
                End If
            End If

            If T > DTL.MathEx.Common.Max(VTc, Vz) Then
                V = 1
                L = 0
            End If

            i = 0
            Do
                If Vz[i] != 0 Then
                    Vy[i] = Vz[i] * Ki[i]/ ((Ki[i]- 1) * V + 1)
                    Vx[i] = Vy[i] / Ki[i]
                Else
                    Vy[i] = 0
                    Vx[i] = 0
                End If
                i += 1
            Loop Until i = n + 1

            i = 0
            soma_x = 0
            soma_y = 0
            Do
                soma_x = soma_x + Vx[i]
                soma_y = soma_y + Vy[i]
                i = i + 1
            Loop Until i = n + 1
            i = 0
            Do
                Vx[i] = Vx[i] / soma_x
                Vy[i] = Vy[i] / soma_y
                i = i + 1
            Loop Until i = n + 1

            Kb_ = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T_, P))
            Kb = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T, P))
            Kb0 = Kb_

            B = Log(Kb_ / Kb) / (1 / T_ - 1 / T)
            A = Log(Kb) - B * (1 / T - 1 / T_)

            for( i = 0;i< n;i++)
                ui[i] = Log(Ki[i]/ Kb)
            Next

            If My.Settings.EnableParallelProcessing Then
                My.MyApplication.IsRunningParallelTasks = True
                Try
                    Dim task1 As Task = New Task(Sub()
                                                     DHv1 = PP.DW_CalcEnthalpyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                                                     DHv2 = PP.DW_CalcEnthalpyDeparture(Vy, Tref, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                                                     C = DHv2
                                                     D = (DHv1 - C) / (T - Tref)
                                                 End Sub)
                    Dim task2 As Task = New Task(Sub()
                                                     If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                                                         DHl1 = PP.DW_CalcEnthalpyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                                                         DHl2 = PP.DW_CalcEnthalpyDeparture(Vx, Tref, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                                                         E = DHl2
                                                         F = (DHl1 - E) / (T - Tref)
                                                     Else
                                                         E = 0
                                                         F = 0
                                                     End If
                                                 End Sub)
                    task1.Start()
                    task2.Start()
                    Task.WaitAll(task1, task2)
                Catch ae As AggregateException
                    For Each ex As Exception In ae.InnerExceptions
                        Throw
                    Next
                End Try
                My.MyApplication.IsRunningParallelTasks = False
            Else
                DHv1 = PP.DW_CalcEnthalpyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                DHv2 = PP.DW_CalcEnthalpyDeparture(Vy, Tref, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                C = DHv2
                D = (DHv1 - C) / (T - Tref)
                If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                    DHl1 = PP.DW_CalcEnthalpyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                    DHl2 = PP.DW_CalcEnthalpyDeparture(Vx, Tref, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                    E = DHl2
                    F = (DHl1 - E) / (T - Tref)
                Else
                    E = 0
                    F = 0
                End If
            End If

            Dim fx(n + 6), x(n + 6), dfdx(n + 6, n + 6), dx(n + 6), xbr(n + 6), fbr(n + 6) As Double
            Dim bo As New BrentOpt.Brent
            Dim bo2 As New BrentOpt.BrentMinimize

            Dim fr As Double

            ecount = 0
            R = Kb * V / (Kb * V + Kb0 * L)

restart:    Do

                //--------------------------------------------------------------
                // STEPS 2, 3, 4, 5, 6, 7 and 8 - Calculate R and Energy Balance
                //--------------------------------------------------------------

                Rant = R
                Tant = T

                Dim dfr, R0, R1 As Double
                Dim icount As Integer = 0

                Do
                    R0 = R
                    If R > 0.99 Then
                        R1 = R - 0.0001
                        fr = Me.EnergyBalance(R0)
                        dfr = (fr - Me.EnergyBalance(R1)) / 0.0001
                    Else
                        R1 = R + 0.0001
                        fr = Me.EnergyBalance(R0)
                        dfr = (fr - Me.EnergyBalance(R1)) / -0.0001
                    End If
                    R0 = R
                    If Abs(fr) < itol Then Exit Do
                    R += -0.3 * fr / dfr
                    If R < 0 Then R = 0
                    If R > 1 Then R = 1
                    icount += 1
                Loop Until icount > maxit_i Or R = 0 Or R = 1 //Or Abs(R - R0) < 0.000001

                Me.EnergyBalance(R)

                //At this point, we have converged T and R for the simplified model. Proceed to step 9.

                //----------------------------------------------------------
                // STEP 9 - Rigorous model Enthalpy and K-values calculation
                //----------------------------------------------------------

                Ki = PP.DW_CalcKvalue(Vx, Vy, T, P)

                If Abs(R - Rant) > 0.01 And Abs(T - Tant) > 0.01 Then
                    for( i = 0;i< n;i++)
                        If Ki[i]!= 0 Then
                            uic[i] = Log(Ki[i]/ Kb)
                        Else
                            uic[i] = ui[i]
                        End If
                    Next
                End If

                Bc = Log(Kb_ / Kb) / (1 / T_ - 1 / T)
                Ac = Log(Kb) - Bc * (1 / T - 1 / T_)

                If My.Settings.EnableParallelProcessing Then
                    My.MyApplication.IsRunningParallelTasks = True
                    Try
                        Dim task1 As Task = New Task(Sub()
                                                         DHv1 = PP.DW_CalcEnthalpyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                                                         DHv2 = PP.DW_CalcEnthalpyDeparture(Vy, T0, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                                                         Cc = DHv2
                                                         Dc = (DHv1 - Cc) / (T - T0)
                                                     End Sub)
                        Dim task2 As Task = New Task(Sub()
                                                         If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                                                             DHl1 = PP.DW_CalcEnthalpyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                                                             DHl2 = PP.DW_CalcEnthalpyDeparture(Vx, T0, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                                                             Ec = DHl2
                                                             Fc = (DHl1 - Ec) / (T - T0)
                                                         Else
                                                             Ec = 0
                                                             Fc = 0
                                                         End If
                                                     End Sub)
                        task1.Start()
                        task2.Start()
                        Task.WaitAll(task1, task2)
                    Catch ae As AggregateException
                        For Each ex As Exception In ae.InnerExceptions
                            Throw
                        Next
                    End Try
                    My.MyApplication.IsRunningParallelTasks = False
                Else
                    DHv1 = PP.DW_CalcEnthalpyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                    DHv2 = PP.DW_CalcEnthalpyDeparture(Vy, T0, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy) / 1000
                    Cc = DHv2
                    Dc = (DHv1 - Cc) / (T - T0)
                    If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                        DHl1 = PP.DW_CalcEnthalpyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                        DHl2 = PP.DW_CalcEnthalpyDeparture(Vx, T0, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx) / 1000
                        Ec = DHl2
                        Fc = (DHl1 - Ec) / (T - T0)
                    Else
                        Ec = 0
                        Fc = 0
                    End If
                End If

                //-------------------------------------------
                // STEP 10 - Update variables using Broyden
                //-------------------------------------------

                for( i = 0;i< n;i++)
                    fx[i] = (ui[i] - uic[i])
                    x(i) = ui[i]
                Next
                fx(n + 1) = (A - Ac)
                fx(n + 2) = (B - Bc)
                fx(n + 3) = (C - Cc)
                fx(n + 4) = (D - Dc)
                fx(n + 5) = (E - Ec)
                fx(n + 6) = (F - Fc)
                x(n + 1) = A
                x(n + 2) = B
                x(n + 3) = C
                x(n + 4) = D
                x(n + 5) = E
                x(n + 6) = F

                If ecount = 0 Then
                    for( i = 0;i< n;i++) + 6
                        For j = 0 To n + 6
                            If i = j Then dfdx(i, j) = 1 Else dfdx(i, j) = 0
                        Next
                    Next
                    Broyden.broydn(n + 6, x, fx, dx, xbr, fbr, dfdx, 0)
                Else
                    Broyden.broydn(n + 6, x, fx, dx, xbr, fbr, dfdx, 1)
                End If

                Dim bo3 As New BrentOpt.BrentMinimize
                bo3.DefineFuncDelegate(AddressOf MinimizeError)
                Dim alpha As Double = 1.0#, err As Double

                tmpdx = dx.Clone
                currx = x.Clone

                ReDim refx(n + 6)

                for( i = 0;i< n;i++)
                    refx[i] = uic[i]
                Next
                refx(n + 1) = Ac
                refx(n + 2) = Bc
                refx(n + 3) = Cc
                refx(n + 4) = Dc
                refx(n + 5) = Ec
                refx(n + 6) = Fc

                err = 0
                If Not PP._ioquick Then err = bo3.brentoptimize(0.2, 2, 0.0001, alpha)

                for( i = 0;i< n;i++)
                    ui[i] = ui[i] + alpha * dx(i)
                Next
                A += alpha * dx(n + 1)
                B += alpha * dx(n + 2)
                C += alpha * dx(n + 3)
                D += alpha * dx(n + 4)
                E += alpha * dx(n + 5)
                F += alpha * dx(n + 6)

                ecount += 1

                If ecount > maxit_e Then Throw New Exception("The flash algorithm reached the maximum number of external iterations.")
                If Double.IsNaN(AbsSum(fx)) Then Throw New Exception(DTL.App.GetLocalString("PropPack_FlashError"))

                Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", T = " & T)
                Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", VF = " & V)
                Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", H error = " & fr)
                Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", Damping Factor = " & alpha)



            Loop Until SumSqr(fx) < etol

            If Abs(fr) > itol Then
                If V <= 0.01 Then
                    //single phase solution found (liquid only). Obtain T using single phase calculation.
                    Dim x1, fx2, dfdx2 As Double
                    ecount = 0
                    If Tref = 0 Then Tref = 298.15
                    x1 = Tref
                    Do
                        fx2 = EnergyBalanceSPL(x1, Nothing)
                        If Math.Abs(fx2) < etol Then Exit Do
                        dfdx2 = (EnergyBalanceSPL(x1 + 1, Nothing) - fx2)
                        x1 = x1 - fx2 / dfdx2
                        ecount += 1
                    Loop Until ecount > maxit_e Or Double.IsNaN(x1)
                    T = x1
                    Vx = Vz
                Else
                    //single phase solution found (vapor only). Obtain T using single phase calculation.
                    Dim x1, fx2, dfdx2 As Double
                    ecount = 0
                    If Tref = 0 Then Tref = 298.15
                    x1 = Tref
                    Do
                        fx2 = EnergyBalanceSPV(x1, Nothing)
                        If Math.Abs(fx2) < etol Then Exit Do
                        dfdx2 = (EnergyBalanceSPV(x1 + 1, Nothing) - fx2)
                        x1 = x1 - fx2 / dfdx2
                        ecount += 1
                    Loop Until ecount > maxit_e Or Double.IsNaN(x1)
                    T = x1
                    Vy = Vz
                End If
                //confirm single-phase solution with a PT Flash.
                Dim res As Object = Me.Flash_PT(Vz, P, T, PP, False, Nothing)
                If Abs(L - res(0)) > 0.0001 And Abs(V - res(1)) > 0.0001 Then
                    //NOT SP solution. go back to 2-phase loop.
                    GoTo restart
                End If
            End If

            d2 = Date.Now

            dt = d2 - d1


            Console.WriteLine("PH Flash [IO]: Converged in " & ecount & " iterations. Time taken: " & dt.TotalMilliseconds & " ms. Error function value: " & AbsSum(fx))

            Return New Object() {L, V, Vx, Vy, T, ecount, Ki, 0.0#, PP.RET_NullVector, 0.0#, PP.RET_NullVector}

        End Function

        Public Overrides Function Flash_PS(ByVal Vz As Double(), ByVal P As Double, ByVal S As Double, ByVal Tref As Double, ByVal PP As PropertyPackages.PropertyPackage, Optional ByVal ReuseKI As Boolean = False, Optional ByVal PrevKi As Double() = Nothing) As Object

            Dim d1, d2 As Date, dt As TimeSpan

            d1 = Date.Now

            maxit_i = CInt(PP.Parameters("PP_PSFMII"))
            maxit_e = CInt(PP.Parameters("PP_PSFMEI"))
            itol = CDbl(PP.Parameters("PP_PSFILT"))
            etol = CDbl(PP.Parameters("PP_PSFELT"))

            n = UBound(Vz)

            proppack = PP
            Sf = S * PP.AUX_MMM(Vz)
            Pf = P

            ReDim Vn(n), Vx(n), Vy(n), Vp(n), ui(n), uic(n), pi(n), Ki(n), fi(n), Vt(n), Vpc(n), VTc(n), Vw(n)

            Vn = PP.RET_VNAMES()
            VTc = PP.RET_VTC
            fi = Vz.Clone

            Tmin = 0
            Tmax = 0
            If Tref = 0 Or Double.IsNaN(Tref) Then
                i = 0
                Tref = 0
                Do
                    Tref += 0.8 * Vz[i] * VTc(i)
                    Tmin += 0.1 * Vz[i] * VTc(i)
                    Tmax += 2.0 * Vz[i] * VTc(i)
                    i += 1
                Loop Until i = n + 1
            Else
                Tmin = Tref - 200
                Tmax = Tref + 200
            End If
            If Tmin < 100 Then Tmin = 100

            //--------------------------------------
            // STEP 1 - Assume u, A, B, C, D, E, F 
            //--------------------------------------

            T = Tref - 0.1
            T_ = Tref - 0.2
            T0 = Tref

            //----------------------------------------
            // STEP 1.1 - Estimate K, Vx, Vy, V and L 
            //----------------------------------------

            //Calculate Ki`s

            Vpc = PP.RET_VPC
            VTc = PP.RET_VTC
            Vw = PP.RET_VW

            If Not ReuseKI Then
                i = 0
                Do
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    Ki[i]= Vpc(i) / P * Math.Exp(5.373 * (1 + Vw(i)) * (1 - VTc(i) / T))
                    i += 1
                Loop Until i = n + 1
            Else
                for( i = 0;i< n;i++)
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    Ki[i]= PrevKi[i]
                Next
            End If

            //Estimate V

            i = 0
            Px = 0
            Do
                Px = Px + (Vz[i] / Vp[i])
                i = i + 1
            Loop Until i = n + 1
            Px = 1 / Px
            Pmin = Px
            i = 0
            Px = 0
            Do
                Px = Px + Vz[i] * Vp[i]
                i = i + 1
            Loop Until i = n + 1

            Pmax = Px
            Pb = Pmax
            Pd = Pmin

            If P <= Pd Then
                //vapor only
                L = 0
                V = 1
            ElseIf P >= Pb Then
                //liquid only
                L = 1
                V = 0
            Else
                //VLE
                V = 1 - (P - Pd) / (Pb - Pd)
                L = 1 - V
            End If

            If n = 0 Then
                If Vp(0) <= P Then
                    L = 1
                    V = 0
                Else
                    L = 0
                    V = 1
                End If
            End If

            If T > DTL.MathEx.Common.Max(proppack.RET_VTC, Vz) Then
                V = 1
                L = 0
            End If

            i = 0
            Do
                If Vz[i] != 0 Then
                    Vy[i] = Vz[i] * Ki[i]/ ((Ki[i]- 1) * V + 1)
                    Vx[i] = Vy[i] / Ki[i]
                Else
                    Vy[i] = 0
                    Vx[i] = 0
                End If
                i += 1
            Loop Until i = n + 1

            i = 0
            soma_x = 0
            soma_y = 0
            Do
                soma_x = soma_x + Vx[i]
                soma_y = soma_y + Vy[i]
                i = i + 1
            Loop Until i = n + 1
            i = 0
            Do
                Vx[i] = Vx[i] / soma_x
                Vy[i] = Vy[i] / soma_y
                i = i + 1
            Loop Until i = n + 1

            Kb_ = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T_, P))
            Kb0 = Kb_

            Kb = CalcKbj1(Ki)

            B = Log(Kb_ / Kb) / (1 / T_ - 1 / T)
            A = Log(Kb) - B * (1 / T - 1 / T_)

            for( i = 0;i< n;i++)
                ui[i] = Log(Ki[i]/ Kb)
            Next

            If My.Settings.EnableParallelProcessing Then
                My.MyApplication.IsRunningParallelTasks = True
                Try
                    Dim task1 As Task = New Task(Sub()
                                                     DSv1 = PP.DW_CalcEntropyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                                                     DSv2 = PP.DW_CalcEntropyDeparture(Vy, Tref, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                                                     C = DSv2
                                                     D = (DSv1 - C) / (T - Tref)
                                                 End Sub)
                    Dim task2 As Task = New Task(Sub()
                                                     If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                                                         DSl1 = PP.DW_CalcEntropyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                                                         DSl2 = PP.DW_CalcEntropyDeparture(Vx, Tref, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                                                         E = DSl2
                                                         F = (DSl1 - E) / (T - Tref)
                                                     Else
                                                         E = 0
                                                         F = 0
                                                     End If
                                                 End Sub)
                    task1.Start()
                    task2.Start()
                    Task.WaitAll(task1, task2)
                Catch ae As AggregateException
                    For Each ex As Exception In ae.InnerExceptions
                        Throw
                    Next
                End Try
                My.MyApplication.IsRunningParallelTasks = False
            Else
                DSv1 = PP.DW_CalcEntropyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                DSv2 = PP.DW_CalcEntropyDeparture(Vy, Tref, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                C = DSv2
                D = (DSv1 - C) / (T - Tref)
                If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                    DSl1 = PP.DW_CalcEntropyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                    DSl2 = PP.DW_CalcEntropyDeparture(Vx, Tref, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                    E = DSl2
                    F = (DSl1 - E) / (T - Tref)
                Else
                    E = 0
                    F = 0
                End If
            End If

            Dim fx(n + 6), x(n + 6), dfdx(n + 6, n + 6), dx(n + 6), xbr(n + 6), fbr(n + 6) As Double
            Dim bo As New BrentOpt.Brent
            Dim bo2 As New BrentOpt.BrentMinimize

            Dim fr As Double

            ecount = 0
            R = Kb * V / (Kb * V + Kb0 * L)

restart:    Do

                //--------------------------------------------------------------
                // STEPS 2, 3, 4, 5, 6, 7 and 8 - Calculate R and Entropy Balance
                //--------------------------------------------------------------

                Rant = R
                Tant = T


                Dim dfr, R0, R1 As Double
                Dim icount As Integer = 0

                Do
                    R0 = R
                    If R > 0.999 Then
                        R1 = R - 0.001
                        fr = Me.EntropyBalance(R0)
                        dfr = (fr - Me.EntropyBalance(R1)) / 0.001
                    Else
                        R1 = R + 0.001
                        fr = Me.EntropyBalance(R0)
                        dfr = (fr - Me.EntropyBalance(R1)) / -0.001
                    End If
                    R0 = R
                    If Abs(fr) < itol Then Exit Do
                    R += -0.3 * fr / dfr
                    If R < 0 Then R = 0
                    If R > 1 Then R = 1
                    icount += 1
                Loop Until icount > maxit_i Or R = 0 Or R = 1 //Or Abs(R - R0) < 0.000001

                Me.EntropyBalance(R)

                //At this point, we have converged T and R for the simplified model. Proceed to step 9.

                //----------------------------------------------------------
                // STEP 9 - Rigorous model Entropy and K-values calculation
                //----------------------------------------------------------

                Ki = PP.DW_CalcKvalue(Vx, Vy, T, P)

                If Abs(R - Rant) > 0.00001 And Abs(T - Tant) > 0.01 Then
                    for( i = 0;i< n;i++)
                        uic[i] = Log(Ki[i]/ Kb)
                    Next
                End If

                Bc = Log(Kb_ / Kb) / (1 / T_ - 1 / T)
                Ac = Log(Kb) - Bc * (1 / T - 1 / T_)

                If My.Settings.EnableParallelProcessing Then
                    My.MyApplication.IsRunningParallelTasks = True
                    Try
                        Dim task1 As Task = New Task(Sub()
                                                         DSv1 = PP.DW_CalcEntropyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                                                         DSv2 = PP.DW_CalcEntropyDeparture(Vy, T0, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                                                         Cc = DSv2
                                                         Dc = (DSv1 - Cc) / (T - T0)
                                                     End Sub)
                        Dim task2 As Task = New Task(Sub()
                                                         If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                                                             DSl1 = PP.DW_CalcEntropyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                                                             DSl2 = PP.DW_CalcEntropyDeparture(Vx, T0, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                                                             Ec = DSl2
                                                             Fc = (DSl1 - Ec) / (T - T0)
                                                         Else
                                                             Ec = 0
                                                             Fc = 0
                                                         End If
                                                     End Sub)
                        task1.Start()
                        task2.Start()
                        Task.WaitAll(task1, task2)
                    Catch ae As AggregateException
                        For Each ex As Exception In ae.InnerExceptions
                            Throw
                        Next
                    End Try
                    My.MyApplication.IsRunningParallelTasks = False
                Else
                    DSv1 = PP.DW_CalcEntropyDeparture(Vy, T, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                    DSv2 = PP.DW_CalcEntropyDeparture(Vy, T0, P, PropertyPackages.State.Vapor) * PP.AUX_MMM(Vy)
                    Cc = DSv2
                    Dc = (DSv1 - Cc) / (T - T0)
                    If T < DTL.MathEx.Common.Max(VTc, Vz) Then
                        DSl1 = PP.DW_CalcEntropyDeparture(Vx, T, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                        DSl2 = PP.DW_CalcEntropyDeparture(Vx, T0, P, PropertyPackages.State.Liquid) * PP.AUX_MMM(Vx)
                        Ec = DSl2
                        Fc = (DSl1 - Ec) / (T - T0)
                    Else
                        Ec = 0
                        Fc = 0
                    End If
                End If

                //-------------------------------------------
                // STEP 10 - Update variables using Broyden
                //-------------------------------------------

                for( i = 0;i< n;i++)
                    fx[i] = (ui[i] - uic[i])
                    x(i) = ui[i]
                Next
                fx(n + 1) = (A - Ac)
                fx(n + 2) = (B - Bc)
                fx(n + 3) = (C - Cc)
                fx(n + 4) = (D - Dc)
                fx(n + 5) = (E - Ec)
                fx(n + 6) = (F - Fc)
                x(n + 1) = A
                x(n + 2) = B
                x(n + 3) = C
                x(n + 4) = D
                x(n + 5) = E
                x(n + 6) = F

                If ecount = 0 Then
                    for( i = 0;i< n;i++) + 6
                        For j = 0 To n + 6
                            If i = j Then dfdx(i, j) = 1 Else dfdx(i, j) = 0
                        Next
                    Next
                    Broyden.broydn(n + 6, x, fx, dx, xbr, fbr, dfdx, 0)
                Else
                    Broyden.broydn(n + 6, x, fx, dx, xbr, fbr, dfdx, 1)
                End If

                Dim bo3 As New BrentOpt.BrentMinimize
                bo3.DefineFuncDelegate(AddressOf MinimizeError)
                Dim alpha As Double = 1.0#, err As Double

                tmpdx = dx.Clone
                currx = x.Clone

                ReDim refx(n + 6)

                for( i = 0;i< n;i++)
                    refx[i] = uic[i]
                Next
                refx(n + 1) = Ac
                refx(n + 2) = Bc
                refx(n + 3) = Cc
                refx(n + 4) = Dc
                refx(n + 5) = Ec
                refx(n + 6) = Fc

                err = 0
                If Not PP._ioquick Then err = bo3.brentoptimize(0, 2, 0.0001, alpha)

                for( i = 0;i< n;i++)
                    ui[i] = ui[i] + alpha * dx(i)
                Next
                A += alpha * dx(n + 1)
                B += alpha * dx(n + 2)
                C += alpha * dx(n + 3)
                D += alpha * dx(n + 4)
                E += alpha * dx(n + 5)
                F += alpha * dx(n + 6)

                ecount += 1

                If ecount > maxit_e Then Throw New Exception("The flash algorithm reached the maximum number of external iterations.")
                If Double.IsNaN(AbsSum(fx)) Then Throw New Exception(DTL.App.GetLocalString("PropPack_FlashError"))

                Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", T = " & T)
                Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", VF = " & V)
                Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", H error = " & fr)
                Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", Damping Factor = " & alpha)



            Loop Until AbsSum(fx) < etol

            If Abs(fr) > itol Then
                If V <= 0.01 Then
                    //single phase solution found (liquid only). Obtain T using single phase calculation.
                    Dim x1, fx2, dfdx2 As Double
                    ecount = 0
                    If Tref = 0 Then Tref = 298.15
                    x1 = Tref
                    Do
                        fx2 = EntropyBalanceSPL(x1, Nothing)
                        If Math.Abs(fx2) < etol Then Exit Do
                        dfdx2 = (EntropyBalanceSPL(x1 + 1, Nothing) - fx2)
                        x1 = x1 - fx2 / dfdx2
                        ecount += 1
                    Loop Until ecount > maxit_e Or Double.IsNaN(x1)
                    T = x1
                    Vx = Vz
                Else
                    //single phase solution found (vapor only). Obtain T using single phase calculation.
                    Dim x1, fx2, dfdx2 As Double
                    ecount = 0
                    If Tref = 0 Then Tref = 298.15
                    x1 = Tref
                    Do
                        fx2 = EntropyBalanceSPV(x1, Nothing)
                        If Math.Abs(fx2) < etol Then Exit Do
                        dfdx2 = (EntropyBalanceSPV(x1 + 1, Nothing) - fx2)
                        x1 = x1 - fx2 / dfdx2
                        ecount += 1
                    Loop Until ecount > maxit_e Or Double.IsNaN(x1)
                    T = x1
                    Vy = Vz
                End If
                //confirm single-phase solution with a PT Flash.
                Dim res As Object = Me.Flash_PT(Vz, P, T, PP, False, Nothing)
                If Abs(L - res(0)) > 0.0001 And Abs(V - res(1)) > 0.0001 Then
                    //NOT single-phase solution. go back to 2-phase loop.
                    GoTo restart
                End If
            End If

            d2 = Date.Now

            dt = d2 - d1

            Console.WriteLine("PS Flash [IO]: Converged in " & ecount & " iterations. Time taken: " & dt.TotalMilliseconds & " ms. Error function value: " & AbsSum(fx))

            Return New Object() {L, V, Vx, Vy, T, ecount, Ki, 0.0#, PP.RET_NullVector, 0.0#, PP.RET_NullVector}

        End Function

        Public Overrides Function Flash_PV(ByVal Vz As Double(), ByVal P As Double, ByVal V As Double, ByVal Tref As Double, ByVal PP As PropertyPackages.PropertyPackage, Optional ByVal ReuseKI As Boolean = False, Optional ByVal PrevKi As Double() = Nothing) As Object

            Dim d1, d2 As Date, dt As TimeSpan

            d1 = Date.Now

            etol = CDbl(PP.Parameters("PP_PTFELT"))
            maxit_e = CInt(PP.Parameters("PP_PTFMEI"))
            itol = CDbl(PP.Parameters("PP_PTFILT"))
            maxit_i = CInt(PP.Parameters("PP_PTFMII"))

            n = UBound(Vz)

            proppack = PP
            Vf = V
            L = 1 - V
            Lf = 1 - Vf
            Pf = P

            ReDim Vn(n), Vx(n), Vy(n), Vp(n), ui(n), uic(n), pi(n), Ki(n), fi(n)
            Dim Vt(n), VTc(n), Tmin, Tmax As Double

            Vn = PP.RET_VNAMES()
            VTc = PP.RET_VTC
            fi = Vz.Clone

            //--------------------------------------
            // STEP 1 - Assume u, A, B, C, D, E, F 
            //--------------------------------------

            Tmin = 0
            Tmax = 0
            If Tref = 0 Then
                i = 0
                Tref = 0
                Do
                    Tref += 0.8 * Vz[i] * VTc(i)
                    Tmin += 0.1 * Vz[i] * VTc(i)
                    Tmax += 2.0 * Vz[i] * VTc(i)
                    i += 1
                Loop Until i = n + 1
            Else
                Tmin = Tref - 50
                Tmax = Tref + 50
            End If

            Dim fx(n + 2), x(n + 2), dfdx(n + 2, n + 2), dx(n + 2), xbr(n + 2), fbr(n + 2) As Double

            T = Tref
            T_ = T - 2
            T0 = T - 5

            If PP.AUX_IS_SINGLECOMP(Vz) Then
                i = 0
                T = 0
                Do
                    T += Vz[i] * PP.AUX_TSATi(P, i)
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    Ki[i]= Vp[i] / P
                    i += 1
                Loop Until i = n + 1
                Vx = Vz
                Vy = Vz
                GoTo final
            End If

            //----------------------------------------
            // STEP 1.1 - Estimate K, Vx, Vy, V and L 
            //----------------------------------------

            //Calculate Ki`s

            If Not ReuseKI Then
                i = 0
                Do
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    Ki[i]= Vp[i] / P
                    i += 1
                Loop Until i = n + 1
            Else
                If Not PP.AUX_CheckTrivial(PrevKi) Then
                    for( i = 0;i< n;i++)
                        Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                        Ki[i]= PrevKi[i]
                    Next
                Else
                    i = 0
                    Do
                        Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                        Ki[i]= Vp[i] / P
                        i += 1
                    Loop Until i = n + 1
                End If
            End If

            i = 0
            Do
                If Vz[i] != 0 Then
                    Vy[i] = Vz[i] * Ki[i]/ ((Ki[i]- 1) * V + 1)
                    Vx[i] = Vy[i] / Ki[i]
                Else
                    Vy[i] = 0
                    Vx[i] = 0
                End If
                i += 1
            Loop Until i = n + 1

            i = 0
            soma_x = 0
            soma_y = 0
            Do
                soma_x = soma_x + Vx[i]
                soma_y = soma_y + Vy[i]
                i = i + 1
            Loop Until i = n + 1
            i = 0
            Do
                Vx[i] = Vx[i] / soma_x
                Vy[i] = Vy[i] / soma_y
                i = i + 1
            Loop Until i = n + 1

            Kb_ = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T_, P))
            Kb0 = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T0, P))
            Kb = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T, P))

            B = Log(Kb_ / Kb) / (1 / T_ - 1 / T)
            A = Log(Kb) - B * (1 / T - 1 / T0)

            for( i = 0;i< n;i++)
                ui[i] = Log(Ki[i]/ Kb)
            Next

            Dim RLoop As Boolean = True

            If V = 0.0# Or V = 1.0# Then RLoop = False

            ecount = 0

            Do

                //--------------------------------------------------------------
                // STEPS 2, 3, 4, 5, 6, 7 and 8 - Calculate R and Energy Balance
                //--------------------------------------------------------------

                R = Kb * V / (Kb * V + Kb0 * L)

                If RLoop Then

                    Dim fr, dfr, R0, R1 As Double
                    Dim icount As Integer = 0

                    Do
                        R1 = R + 0.001
                        fr = Me.LiquidFractionBalance(R)
                        dfr = (fr - Me.LiquidFractionBalance(R1)) / -0.001
                        R0 = R
                        R += -fr / dfr
                        If R < 0 Then R = 0
                        If R > 1 Then R = 1
                        icount += 1
                    Loop Until Abs(fr) < itol Or icount > maxit_i Or Abs(R - R0) < 0.000001

                Else

                    for( i = 0;i< n;i++)
                        pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
                    Next

                    Dim sumpi As Double = 0
                    Dim sumeuipi As Double = 0

                    for( i = 0;i< n;i++)
                        sumpi += pi[i]
                        sumeuipi += Math.Exp(ui[i]) * pi[i]
                    Next

                    Kb = sumpi / sumeuipi

                    Dim Tant As Double = T
                    T = 1 / T0 + (Log(Kb) - A) / B
                    T = 1 / T

                    for( i = 0;i< n;i++)
                        Vx[i] = pi[i] / sumpi
                        Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
                    Next

                End If

                //At this point, we have converged T for the simplified model. Proceed to step 9.

                //----------------------------------------------------------
                // STEP 9 - Rigorous model Enthalpy and K-values calculation
                //----------------------------------------------------------

                Ki = PP.DW_CalcKvalue(Vx, Vy, T, P)
                //Kb_ = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T_, P))

                for( i = 0;i< n;i++)
                    uic[i] = Log(Ki[i]/ Kb)
                Next

                Bc = Log(Kb_ / Kb) / (1 / T_ - 1 / T)
                Ac = Log(Kb) - Bc * (1 / T - 1 / T0)

                //-------------------------------------------
                // STEP 10 - Update variables using Broyden
                //-------------------------------------------

                for( i = 0;i< n;i++)
                    fx[i] = (ui[i] - uic[i])
                    x(i) = ui[i]
                Next
                fx(n + 1) = (A - Ac)
                fx(n + 2) = (B - Bc)
                x(n + 1) = A
                x(n + 2) = B

                If ecount = 0 Then
                    for( i = 0;i< n;i++) + 2
                        For j = 0 To n + 2
                            If i = j Then dfdx(i, j) = 1 Else dfdx(i, j) = 0
                        Next
                    Next
                    Broyden.broydn(n + 2, x, fx, dx, xbr, fbr, dfdx, 0)
                Else
                    Broyden.broydn(n + 2, x, fx, dx, xbr, fbr, dfdx, 1)
                End If


                Dim bo2 As New BrentOpt.BrentMinimize
                bo2.DefineFuncDelegate(AddressOf MinimizeError)
                Dim alpha As Double = 1.0#, err As Double

                tmpdx = dx.Clone
                currx = x.Clone

                ReDim refx(n + 2)

                for( i = 0;i< n;i++)
                    refx[i] = uic[i]
                Next
                refx(n + 1) = Ac
                refx(n + 2) = Bc

                err = 0
                If Not PP._ioquick Then err = bo2.brentoptimize(0, 2, 0.0001, alpha)

                for( i = 0;i< n;i++)
                    ui[i] = ui[i] + alpha * dx(i)
                Next
                A += alpha * dx(n + 1)
                B += alpha * dx(n + 2)

                //for( i = 0;i< n;i++)
                //    ui[i] = uic[i]
                //Next
                //A = Ac
                //B = Bc

                ecount += 1

                If ecount > maxit_e Then
                    Throw New Exception("The flash algorithm reached the maximum number of external iterations.")
                End If
                If Double.IsNaN(AbsSum(fx)) Then
                    Throw New Exception(DTL.App.GetLocalString("PropPack_FlashError"))
                End If

                Console.WriteLine("PV Flash [IO]: Iteration #" & ecount & ", T = " & T & ", VF = " & V)
                Console.WriteLine("PV Flash [IO]: Iteration #" & ecount & ", Damping Factor = " & alpha)



            Loop Until AbsSum(fx) < etol * (n + 2)

final:      d2 = Date.Now

            dt = d2 - d1

            Console.WriteLine("PV Flash [IO]: Converged in " & ecount & " iterations. Time taken: " & dt.TotalMilliseconds & " ms. Error function value: " & AbsSum(fx))

            Return New Object() {L, V, Vx, Vy, T, ecount, Ki, 0.0#, PP.RET_NullVector, 0.0#, PP.RET_NullVector}

        End Function

        Public Overrides Function Flash_TV(ByVal Vz As Double(), ByVal T As Double, ByVal V As Double, ByVal Pref As Double, ByVal PP As PropertyPackages.PropertyPackage, Optional ByVal ReuseKI As Boolean = False, Optional ByVal PrevKi As Double() = Nothing) As Object

            Dim d1, d2 As Date, dt As TimeSpan

            d1 = Date.Now

            etol = CDbl(PP.Parameters("PP_PTFELT"))
            maxit_e = CInt(PP.Parameters("PP_PTFMEI"))

            n = UBound(Vz)

            proppack = PP
            Vf = V
            L = 1 - V
            Lf = 1 - Vf
            Tf = T

            ReDim Vn(n), Vx(n), Vy(n), Vp(n), ui(n), uic(n), pi(n), Ki(n), fi(n)

            Dim VTc = PP.RET_VTC()

            Vn = PP.RET_VNAMES()
            fi = Vz.Clone

            If Pref = 0 Then
                i = 0
                Do
                    If T / VTc(i) <= 0.9 Then
                        Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    End If
                    i += 1
                Loop Until i = n + 1
                Pmin = Common.Min(Vp)
                Pmax = Common.Max(Vp)
                Pref = Pmin + (1 - V) * (Pmax - Pmin)
            Else
                Pmin = Pref * 0.8
                Pmax = Pref * 1.2
            End If

            Dim fx(n + 2), x(n + 2), dfdx(n + 2, n + 2), dx(n + 2), xbr(n + 2), fbr(n + 2) As Double

            P = Pref
            P_ = Pref * 1.05
            P0 = Pref * 0.95

            If PP.AUX_IS_SINGLECOMP(Vz) Then
                i = 0
                P = 0
                Do
                    P += Vz[i] * Vp[i]
                    i += 1
                Loop Until i = n + 1
                Vx = Vz
                Vy = Vz
                GoTo final
            End If

            //----------------------------------------
            // STEP 1.1 - Estimate K, Vx, Vy, V and L 
            //----------------------------------------

            //Calculate Ki`s

            If Not ReuseKI Then
                i = 0
                Do
                    Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                    Ki[i]= Vp[i] / P
                    i += 1
                Loop Until i = n + 1
            Else
                If Not PP.AUX_CheckTrivial(PrevKi) And Not Double.IsNaN(PrevKi(0)) Then
                    for( i = 0;i< n;i++)
                        Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                        Ki[i]= PrevKi[i]
                    Next
                Else
                    i = 0
                    Do
                        Vp[i] = PP.AUX_PVAPi(Vn(i), T)
                        Ki[i]= Vp[i] / P
                        i += 1
                    Loop Until i = n + 1
                End If
            End If

            i = 0
            Do
                If Vz[i] != 0 Then
                    Vy[i] = Vz[i] * Ki[i]/ ((Ki[i]- 1) * V + 1)
                    Vx[i] = Vy[i] / Ki[i]
                Else
                    Vy[i] = 0
                    Vx[i] = 0
                End If
                i += 1
            Loop Until i = n + 1

            i = 0
            soma_x = 0
            soma_y = 0
            Do
                soma_x = soma_x + Vx[i]
                soma_y = soma_y + Vy[i]
                i = i + 1
            Loop Until i = n + 1
            i = 0
            Do
                Vx[i] = Vx[i] / soma_x
                Vy[i] = Vy[i] / soma_y
                i = i + 1
            Loop Until i = n + 1

            Kb0 = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T, P0))
            Kb_ = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T, P_))
            Kb = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T, P))

            B = Log(Kb_ * P_ / (Kb0 * P0)) / Log(P_ / P0)
            A = Log(Kb * P) - B * Log(P / P0)

            for( i = 0;i< n;i++)
                ui[i] = Log(Ki[i]/ Kb)
            Next

            Dim RLoop As Boolean = True

            If V = 0.0# Or V = 1.0# Then RLoop = False

            ecount = 0

            Do

                //--------------------------------------------------------------
                // STEPS 2, 3, 4, 5, 6, 7 and 8 - Calculate R and Energy Balance
                //--------------------------------------------------------------

                R = Kb * V / (Kb * V + Kb0 * L)

                If RLoop Then

                    Dim fr, dfr, R0, R1 As Double
                    Dim icount As Integer = 0

                    Do
                        R1 = R + 0.001
                        fr = Me.LiquidFractionBalanceP(R)
                        dfr = (fr - Me.LiquidFractionBalanceP(R1)) / -0.001
                        R0 = R
                        R += -fr / dfr
                        If R < 0 Then R = 0
                        If R > 1 Then R = 1
                        icount += 1
                    Loop Until Abs(fr) < itol Or icount > maxit_i

                Else

                    for( i = 0;i< n;i++)
                        pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
                    Next

                    Dim sumpi As Double = 0
                    Dim sumeuipi As Double = 0

                    for( i = 0;i< n;i++)
                        sumpi += pi[i]
                        sumeuipi += Math.Exp(ui[i]) * pi[i]
                    Next

                    Kb = sumpi / sumeuipi

                    P = Math.Exp((A - Log(Kb) - B * Log(P0)) / (1 - B))

                    //If Double.IsNaN(P) Then P = Pref * 1.01
                    //If P < Pmin * 0.8 Then
                    //    P = Pmin * 0.8 + 101325
                    //ElseIf P > Pmax * 1.2 Then
                    //    P = Pmax * 1.2 - 101325
                    //End If

                    for( i = 0;i< n;i++)
                        Vx[i] = pi[i] / sumpi
                        Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
                    Next

                End If

                //At this point, we have converged T for the simplified model. Proceed to step 9.

                //----------------------------------------------------------
                // STEP 9 - Rigorous model Enthalpy and K-values calculation
                //----------------------------------------------------------

                Ki = PP.DW_CalcKvalue(Vx, Vy, T, P)
                //Kb_ = CalcKbj1(PP.DW_CalcKvalue(Vx, Vy, T, P_))

                for( i = 0;i< n;i++)
                    uic[i] = Log(Ki[i]/ Kb)
                Next

                Bc = Log(Kb_ * P_ / (Kb0 * P0)) / Log(P_ / P0)
                Ac = Log(Kb * P) - Bc * Log(P / P0)

                //-------------------------------------------
                // STEP 10 - Update variables using Broyden
                //-------------------------------------------

                for( i = 0;i< n;i++)
                    fx[i] = (ui[i] - uic[i])
                    x(i) = ui[i]
                Next
                fx(n + 1) = (A - Ac)
                fx(n + 2) = (B - Bc)
                x(n + 1) = A
                x(n + 2) = B

                If ecount = 0 Then
                    for( i = 0;i< n;i++) + 2
                        For j = 0 To n + 2
                            If i = j Then dfdx(i, j) = 1 Else dfdx(i, j) = 0
                        Next
                    Next
                    Broyden.broydn(n + 2, x, fx, dx, xbr, fbr, dfdx, 0)
                Else
                    Broyden.broydn(n + 2, x, fx, dx, xbr, fbr, dfdx, 1)
                End If

                Dim bo2 As New BrentOpt.BrentMinimize
                bo2.DefineFuncDelegate(AddressOf MinimizeError)
                Dim alpha As Double = 1.0#, err As Double

                tmpdx = dx.Clone
                currx = x.Clone

                ReDim refx(n + 2)

                for( i = 0;i< n;i++)
                    refx[i] = uic[i]
                Next
                refx(n + 1) = Ac
                refx(n + 2) = Bc

                err = 0
                If Not PP._ioquick Then err = bo2.brentoptimize(0, 2, 0.0001, alpha)

                for( i = 0;i< n;i++)
                    ui[i] = ui[i] + alpha * dx(i)
                Next
                A += alpha * dx(n + 1)
                B += alpha * dx(n + 2)

                ecount += 1

                If ecount > maxit_e Then Throw New Exception("The flash algorithm reached the maximum number of external iterations.")
                If Double.IsNaN(AbsSum(fx)) Then Throw New Exception(DTL.App.GetLocalString("PropPack_FlashError"))

                Console.WriteLine("TV Flash [IO]: Iteration #" & ecount & ", P = " & P & ", VF = " & V)
                Console.WriteLine("TV Flash [IO]: Iteration #" & ecount & ", Damping Factor = " & alpha)



            Loop Until AbsSum(fx) < etol * (n + 2)

final:      d2 = Date.Now

            dt = d2 - d1

            Console.WriteLine("TV Flash [IO]: Converged in " & ecount & " iterations. Time taken: " & dt.TotalMilliseconds & " ms. Error function value: " & AbsSum(fx))

            Return New Object() {L, V, Vx, Vy, P, ecount, Ki, 0.0#, PP.RET_NullVector, 0.0#, PP.RET_NullVector}

        End Function

        Private Function LiquidFractionBalance(ByVal R As Double) As Double

            for( i = 0;i< n;i++)
                pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
            Next

            Dim sumpi As Double = 0
            Dim sumeuipi As Double = 0

            for( i = 0;i< n;i++)
                sumpi += pi[i]
                sumeuipi += Math.Exp(ui[i]) * pi[i]
            Next

            Kb = sumpi / sumeuipi

            Dim Tant As Double = T
            T = 1 / T0 + (Log(Kb) - A) / B
            T = 1 / T

            //If T < Tmin Then T = Tmin
            //If T > Tmax Then T = Tmax

            for( i = 0;i< n;i++)
                Vx[i] = pi[i] / sumpi
                Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
            Next

            L = (1 - R) * sumpi
            V = 1 - L

            Dim eberror As Double = L - Lf



            Return eberror

        End Function

        Private Function LiquidFractionBalanceP(ByVal R As Double) As Double

            for( i = 0;i< n;i++)
                pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
            Next

            Dim sumpi As Double = 0
            Dim sumeuipi As Double = 0

            for( i = 0;i< n;i++)
                sumpi += pi[i]
                sumeuipi += Math.Exp(ui[i]) * pi[i]
            Next

            Kb = sumpi / sumeuipi

            P = Math.Exp((A - Log(Kb) - B * Log(P0)) / (1 - B))

            for( i = 0;i< n;i++)
                Vx[i] = pi[i] / sumpi
                Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
            Next

            L = (1 - R) * sumpi
            V = 1 - L

            Dim eberror As Double = L - Lf



            Return eberror

        End Function

        Private Function EnergyBalance(ByVal R As Double) As Double

            for( i = 0;i< n;i++)
                pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
            Next

            Dim sumpi As Double = 0
            Dim sumeuipi As Double = 0

            for( i = 0;i< n;i++)
                sumpi += pi[i]
                sumeuipi += Math.Exp(ui[i]) * pi[i]
            Next

            Kb = sumpi / sumeuipi

            T = 1 / T_ + (Log(Kb) - A) / B
            T = 1 / T

            If T < Tmin Then T = Tmin
            If T > Tmax Then T = Tmax

            for( i = 0;i< n;i++)
                Vx[i] = pi[i] / sumpi
                Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
            Next

            L = (1 - R) * sumpi
            V = 1 - L

            DHv = C + D * (T - T0)
            DHl = E + F * (T - T0)

            If Double.IsNaN(DHl) Then DHl = 0

            Hvid = proppack.RET_Hid(298.15, T, Vy) * proppack.AUX_MMM(Vy) / 1000
            Hlid = proppack.RET_Hid(298.15, T, Vx) * proppack.AUX_MMM(Vx) / 1000

            Dim eberror As Double = Hf / 1000 - L * (DHl + Hlid) - V * (DHv + Hvid)



            Return eberror

        End Function

        Private Function EnergyBalanceAbs(ByVal R As Double) As Double

            for( i = 0;i< n;i++)
                pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
            Next

            Dim sumpi As Double = 0
            Dim sumeuipi As Double = 0

            for( i = 0;i< n;i++)
                sumpi += pi[i]
                sumeuipi += Math.Exp(ui[i]) * pi[i]
            Next

            Kb = sumpi / sumeuipi

            T = 1 / T_ + (Log(Kb) - A) / B
            T = 1 / T

            If T < Tmin Then T = Tmin
            If T > Tmax Then T = Tmax

            for( i = 0;i< n;i++)
                Vx[i] = pi[i] / sumpi
                Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
            Next

            L = (1 - R) * sumpi
            V = 1 - L

            DHv = C + D * (T - T0)
            DHl = E + F * (T - T0)

            Dim sumzihv0i As Double = 0

            for( i = 0;i< n;i++)
                sumzihv0i += fi[i] * proppack.AUX_INT_CPDTi(298.15, T, Vn(i)) * proppack.AUX_MMM(fi) / 1000
            Next

            Dim eberror As Double = L * (DHv - DHl) - sumzihv0i - DHv + Hf / 1000


            Return Abs(eberror)

        End Function

        Private Function EnergyBalanceSPL(ByVal T As Double, ByVal otherargs As Object) As Double



            Dim HL, balerror As Double

            HL = proppack.DW_CalcEnthalpy(fi, T, Pf, PropertyPackages.State.Liquid) * proppack.AUX_MMM(fi)

            balerror = Hf - HL

            Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", T = " & T)
            Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", VF = 0 (SP)")
            Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", H error = " & balerror)

            Return balerror


        End Function

        Private Function EnergyBalanceSPV(ByVal T As Double, ByVal otherargs As Object) As Double



            Dim HV, balerror As Double

            HV = proppack.DW_CalcEnthalpy(fi, T, Pf, PropertyPackages.State.Vapor) * proppack.AUX_MMM(fi)

            balerror = Hf - HV

            Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", T = " & T)
            Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", VF = 1 (SP)")
            Console.WriteLine("PH Flash [IO]: Iteration #" & ecount & ", H error = " & balerror)

            Return balerror


        End Function

        Private Function EntropyBalance(ByVal R As Double) As Double

            for( i = 0;i< n;i++)
                pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
            Next

            Dim sumpi As Double = 0
            Dim sumeuipi As Double = 0

            for( i = 0;i< n;i++)
                sumpi += pi[i]
                sumeuipi += Math.Exp(ui[i]) * pi[i]
            Next

            Kb = sumpi / sumeuipi

            T = 1 / T_ + (Log(Kb) - A) / B
            T = 1 / T

            If T < Tmin Then T = Tmin
            If T > Tmax Then T = Tmax

            for( i = 0;i< n;i++)
                Vx[i] = pi[i] / sumpi
                Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
            Next

            L = (1 - R) * sumpi
            V = 1 - L

            DSv = C + D * (T - T0)
            DSl = E + F * (T - T0)

            Svid = proppack.RET_Sid(298.15, T, Pf, Vy) * proppack.AUX_MMM(Vy)
            Slid = proppack.RET_Sid(298.15, T, Pf, Vx) * proppack.AUX_MMM(Vx)

            Dim eberror As Double = Sf - L * (DSl + Slid) - V * (DSv + Svid)


            Return eberror

        End Function

        Private Function EntropyBalanceAbs(ByVal R As Double) As Double

            for( i = 0;i< n;i++)
                pi[i] = fi[i] / (1 - R + Kb0 * R * Math.Exp(ui[i]))
            Next

            Dim sumpi As Double = 0
            Dim sumeuipi As Double = 0

            for( i = 0;i< n;i++)
                sumpi += pi[i]
                sumeuipi += Math.Exp(ui[i]) * pi[i]
            Next

            Kb = sumpi / sumeuipi

            T = 1 / T_ + (Log(Kb) - A) / B
            T = 1 / T

            If T < Tmin Then T = Tmin
            If T > Tmax Then T = Tmax

            for( i = 0;i< n;i++)
                Vx[i] = pi[i] / sumpi
                Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
            Next

            L = (1 - R) * sumpi
            V = 1 - L

            DSv = C + D * (T - T0)
            DSl = E + F * (T - T0)

            Svid = proppack.RET_Sid(298.15, T, Pf, Vy) * proppack.AUX_MMM(Vy)
            Slid = proppack.RET_Sid(298.15, T, Pf, Vx) * proppack.AUX_MMM(Vx)

            Dim eberror As Double = Sf - L * (DSl + Slid) - V * (DSv + Svid)


            Return Abs(eberror)

        End Function

        Private Function EntropyBalanceSPL(ByVal T As Double, ByVal otherargs As Object) As Double



            Dim SL, balerror As Double

            SL = proppack.DW_CalcEntropy(fi, T, Pf, PropertyPackages.State.Liquid) * proppack.AUX_MMM(fi)

            balerror = Sf - SL

            Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", T = " & T)
            Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", VF = 0 (SP)")
            Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", S error = " & balerror)

            Return balerror


        End Function

        Private Function EntropyBalanceSPV(ByVal T As Double, ByVal otherargs As Object) As Double



            Dim SV, balerror As Double

            SV = proppack.DW_CalcEntropy(fi, T, Pf, PropertyPackages.State.Vapor) * proppack.AUX_MMM(fi)

            balerror = Sf - SV

            Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", T = " & T)
            Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", VF = 1 (SP)")
            Console.WriteLine("PS Flash [IO]: Iteration #" & ecount & ", S error = " & balerror)

            Return balerror

        End Function

        Private Function CalcKbj1(ByVal K() As Double) As Double

            Dim i As Integer
            Dim n As Integer = UBound(K) - 1

            Dim Kbj1 As Object

            Kbj1 = K(0)
            For i = 1 To n
                If Abs(K(i) - 1) < Abs(Kbj1 - 1) Then Kbj1 = K(i)
            Next

            Return Kbj1

        End Function

        Private Function CalcKbj2(ByVal K() As Double, ByVal K2() As Double, ByVal T1 As Double, ByVal T2 As Double) As Double

            Dim i As Integer
            Dim n As Integer = UBound(K) - 1

            Dim Kbj1 As Object
            Dim Kw11(n), Kw21(n) As Object
            Dim wi(n), ti(n), sumti As Double

            for( i = 0;i< n;i++)
                Kw11(i) = K(i)

                Kw21(i) = K2(i)
            Next

            sumti = 0
            for( i = 0;i< n;i++)
                ti(i) = Vy[i] * (Log(Kw21(i)) - Log(Kw11(i))) / (1 / T2 - 1 / T1)
                sumti += ti(i)
            Next

            for( i = 0;i< n;i++)
                wi(i) = ti(i) / sumti
            Next

            Kbj1 = 0
            for( i = 0;i< n;i++)
                Kbj1 += wi(i) * Log(K(i))
                Kbj1 = Math.Exp(Kbj1)
                If Kbj1 < 0 Then
                    Kbj1 = K(i)
                    For j = 1 To n
                        If Abs(K(j) - 1) < Abs(Kbj1 - 1) Then Kbj1 = K(j)
                    Next
                End If
            Next

            Return Kbj1

        End Function

        Private Function TPErrorFunc(ByVal Rt As Double) As Double

            for( i = 0;i< n;i++)
                pi[i] = fi[i] / (1 - Rt + Kb0 * Rt * Math.Exp(ui[i]))
            Next

            Dim sumpi As Double = 0
            Dim sumeuipi As Double = 0

            for( i = 0;i< n;i++)
                sumpi += pi[i]
                sumeuipi += Math.Exp(ui[i]) * pi[i]
            Next

            for( i = 0;i< n;i++)
                Vx[i] = pi[i] / sumpi
                Vy[i] = Math.Exp(ui[i]) * pi[i] / sumeuipi
            Next

            L = (1 - Rt) * sumpi
            V = 1 - L

            Dim eberror As Double = sumpi / sumeuipi - 1



            Return eberror

        End Function

        Private Function MinimizeError(ByVal alpha As Double) As Double

            Dim n As Integer = UBound(tmpdx)
            Dim i As Integer
            Dim errors(n) As Double

            for( i = 0;i< n;i++)
                errors(i) = (refx[i] - (currx(i) + alpha * tmpdx(i))) ^ 2
            Next

            Return Common.Sum(errors)

        End Function

    End Class

End Namespace